Compact, orthogonal, and complete basis sets for solving the Schrodinger equation 
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We present a new type of basis set which is local, compact, and orthogonal. The basis functions, 
called orthlets, are centered at the sites of a lattice and are specifically adapted to represent the 
system being studied. The adaptability includes the ability to have singular behavior within an 
orthlet, allowing a single orthlet to represent a function in the vicinity of a singularity. 



PACS Numbers: 

Most modern numerical solution techniques to the 
Schrodinger equation begin with the introduction of a 
basis set, thereby making an infinite Hilbert space finite. 
Because there are a number of incompatible qualities one 
would like in a basis set, a wide variety of basis sets are 
in use, each of which is better behaved according to some 
set of criteria. Among the desirable qualities are orthog- 
onality, locality, compactness (i.e. compact support), the 
ability to represent space uniformly, the ability to repre- 
sent singular regions with higher resolution, the ability to 
incorporate prior knowledge about singular regions, the 
ability to ignore empty regions, and the availability of 
specialized efficient algorithms (such as the fast fourier 
transforms (FFTs) or wavelet transforms) for doing inte- 
grals and solving differential equations. 

For example, in electronic structure calculations for 
solids using density functional methods, plane waves are 
widely and successfully used [jjj. These are orthogo- 
nal, have uniform resolution, and the FFT allows rapid 
switching between real and fourier space. Pseudopoten- 
tials are normally used to represent atomic cores. How- 
ever, for more accurate treatment of the cores, for non- 
periodic systems (including molecules and surfaces), and 
for more accurate treatment of interatomic correlations, 
the plane wave basis is inconvenient. 

In quantum chemistry, the standard choice for basis 
functions is the product of a radial function centered on 
an atom times a cartesian or spherical harmonic 0] . Be- 
cause the radial functions which solve the Hartree Fock 
equations for atoms are known, remarkably small num- 
bers of basis functions are needed — often only about 
twice as many basis functions as there are electrons. 
The nonorthogonality of the basis is easily dealt with 
in Hartree Fock. The major drawbacks relate to scal- 
ing to large systems and to high accuracy. The number 
of two-electron integrals needed to represent the inter- 
electron Coulomb interaction scales as TV 4 , where N is 
the number of basis functions. Moreover, the orthog- 
onalization required for most treatments of correlations 
beyond Hartree Fock destroys the approximate locality of 
the functions; consequently, computation time typically 
scales as N 6 or worse. 

Wavelet bases are another potentially attractive alter- 
native |3|. These nonorthogonal bases allow for widely 



varying resolution to represent both cores and valence 
electrons. However, hundreds or thousands of wavelets 
on various length scales may be needed to represent a 
second row atomic core, compared to perhaps a dozen of 
the radial basis functions used in quantum chemistry, 

In this letter, we propose a new type of basis set which 
is orthogonal, very localized and compact, which allows 
variable resolution, and which allows prior knowledge 
about singularities to be incorporated into the basis, 
while keeping the number of basis functions to a mini- 
mum. Our approach is most closely related to the finite 
element basis using orthogonal shape functions developed 
by White, Wilkins, and Teter (WWT) f§. The major 
problem with the approach of WWT was the difficulty 
in obtaining adequate resolution for the cores. Our new 
approach overcomes that difficulty. Although these bases 
were developed with electronic structure calculations in 
mind, we expect them to be useful in a variety of other 
contexts as well. 

Consider a set of localized shape functions <fii(r), and 
a lattice {Rj}- We can generate a set of functions for 
each lattice site by translation, <f>ij = <j>i(r — Rj). Let 
the functions have the following properties: 1) the set of 
functions 4>i(r) is orthonormal; 2) each function is also 
orthogonal to the functions on all other lattice sites; and 
3) the total set of functions on all lattice sites is complete. 
(Wannier functions also have these properties.) We de- 
fine a projection operator for site j by 



(i) 



In coordinate notation this operator is 

Pj(r, ?) = E M?- - Rj), (2) 

i 

where the 4>i are assumed real. Completeness implies 
that Ylj Pj = 1- Now consider the application of Pj on 
an arbitrary function f(f): fj{r) — Pjf(r). Then 



and 



(3) 



(4) 
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We will call functions such as fj(r), which are local, or- 
thogonal, and specifically adapted to a function or to a 
set of functions, orthlets. The basis formed by the orth- 
lets fj(r) is in a sense a perfect basis for representing 
f(r): it is orthogonal, represents the function exactly, 
and has the minimum number of functions given the scale 
set by the lattice spacing. We describe below how to form 
an orthlet basis describing an arbitrary set of functions 
f a (r) to a specified accuracy. 

In order for the orthlets to be useful, the shape func- 
tions <fii(r) should be smooth and local, and preferably 
compact. In two or more dimensions, shape functions 
can be written as products of one dimensional shape 
functions [f|, so that we need only consider the ID case. 
WWT developed a set of four shape functions with con- 
tinuous derivatives up to third order, which were able 
to represent exactly polynomials up to third order 
These shape functions were compact, with a total width 
of two, where we assume a lattice spacing of unity hence- 
forth. The compactness means that orthogonality must 
only be specifically arranged for nearest neighbor func- 
tions. Here we give six orthogonal shape functions which 
are smooth, i.e. all derivatives are continuous, which also 
are compact with width two, and which are able to rep- 
resent polynomials exactly up to order five. We believe 
these shape functions are sufficiently complete for most 
uses, because we give alternative ways of generating the 
orthlets in the vicinity of singularities and also for chang- 
ing lattice spacings. The smoothness allows additional 
functions to be added to increase completeness without 
adjusting the functions one already has. 

All shape functions S n (x) are defined for x > 0; S n (x) 
is even (odd) if n is. We construct Sq(x) to represent 
a constant function exactly. We first define a smooth 
"splicing function" p{x) which divides a function to be 
fit to into pieces. We require that p{x) — for |je| > 1, 
p(x) = p{—x), and that 



where 



p(x) +p(x - 1) = 1 

We choose the function 
f 



< x < 1. 



(5) 



p(*) = Htanh(| 0< * <L (6) 

This bell-shaped function has essential singularities at 
x = 0, ±1, allowing it to have compact support yet be 
smooth. 

The shape function is obtained by multiplying the 
function to fit to, in this case unity, by p(x), and then 
adding a smooth oscillating function o(x) to induce or- 
thogonality. The fit will not be spoiled if 



o[x) + o(x - 1) = 0<x<l. 



We choose 



°( x ) — a 0wi°m( a 



(7) 



(8) 



o m {x) = p(2x — 1) sin(2m7r(x — 1/2)). 



(9) 



For the first shape function, Sq(x), is suffices to take 
M = 1, with aoi = -0.507021142747521: 



Sq(x) = p(x) + ao%o\(x) < x < 1. 



(10) 



Note that S (x) + S (x - 1) = 1 for < x < 1, so that 
Sq(x) is already normalized. We show So(x) in Fig. 1(b). 

The second shape function is obtained similarly, by 
requiring an exact fit to the function x. First, we attempt 
to fit the function x using only Sq(x). We then make 
S±(x) out of the error or residual of this fit, 



r\(x) = x — Sq(x — 1) < x < 1. 



(11) 



In this case, the extra orthogonalizing functions must be 
made out of cosine functions rather than sine functions, 
because the function itself is odd. For < x < 1 



Si(x) = Ni ^(x)ri(a;) + ^ o,i m e m (x)j , 

with Si(— x) = — Si(x), 

e m {x) = p(2x - 1) cos((2to - l)n(x - 1/2)), 



(12) 



(13) 



ail = 0.132403793351197, a u = -0.048844623781880, 
and Ni = 3.78750743638139. 
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FIG. 1. (a) and (b) Shape functions S n (x) 
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To determine the coefficients a\ mi and later coefficients 
for higher order shape functions, two different linear com- 
binations of the e m (x) or o m (x) were found which in- 
duced orthogonality to all the lower order shape functions 
centered at x = 1 . These two different combinations were 
then combined in order to induce orthogonality of S n (x) 
and S n (x — 1). This last procedure required solving a 
quadratic equation, which might not have real roots. If 
there were no real roots, more orthogonalizing functions 
were included. A solution in this case required finding 
both positive and negative eigenvalues of a symmetric 
matrix, which also sometimes did not occur. There does 
not appear to be any guarantee of a real solution; in fact, 
we did not find a satisfactory Sq (x) able to fit x 6 exactly. 

The next shape function is given by 

S 2 (x) = N 2 [p(x)(x 2 - c 20 S (x) - d 20 S a (x - 1) (14) 
-d 2 iSi(x - 1)) + a 21 oi(x) + a 22 o 2 (x) + a 23 o 3 (x)}, 

c 20 = 0.0697096675548214, d 20 = 1.0697096675548214, 
d 21 = 0.528051768503122, a 21 = 0.088401702549656, 
a 22 = -0.126644764032427, a 23 = -0.025357986069321, 
and N 2 = 11.9312518524753. Expressions for S 3 to S 5 
will be presented elsewhere ||. 

In Fig. 2, we use these shape functions to generate 
an orthlet basis for a function with a slope discontinu- 
ity. For the lattice sites away from the singularity, the 
overlap integral of the function with each shape function 
was computed numerically. Resumming the shape func- 
tions on a site with these integrals as coefficients gives 
the orthlet. In the case of the orthlet at x = 0, where the 
slope discontinuity is, an expansion would converge too 
slowly. Instead, the orthlet was obtained by subtracting 
from f(x) the orthlets on the adjacent sites. 
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FIG. 2. Orthlets used to fit the function 
f(x) = exp(-|a;|) + l/((x - 3) 2 + 1/2). The sum of these 
four orthlets is equal to f(x) within the region — 1 < x < 2. 



If the shape functions describe the function perfectly in 
the region away from the singularity, then the subtracted 
function will be identically zero for |a;| > 1. If the fit is 
not perfect but very good, the function can be set to zero 
for \x\ > 1. However, a small discontinuity and a small 
lack of orthogonality can result from this procedure. In 
this example, the discontinuities were less than 10~ 5 and 
were ignored, along with a small nonorthogonality. In 
order to ensure perfect continuity and orthogonality in 
this example, we could multiply the subtracted function 
by a smooth windowing function which is unity for most 
of the interval — 1 < x < 1, and is zero for |x| > 1. Then 
the resulting function can be explicitly orthogonalized 
in a Gram-Schmidt (GS) procedure to the neighboring 
shape functions. Because of the orthogonalization, the 
resulting function would extend from — 2 < x < 2. 

Another procedure for dealing with singularities is to 
add a set of localized functions near the singularity, which 
are chosen for their ability to represent the singularity. 
Then an explicit GS orthogonalization procedure mixes 
these functions with the shape functions that overlap 
them. This GS procedure does not destroy overall lo- 
cality: all shape functions not directly overlapping the 
functions added would still automatically be orthogonal 
to the new set of functions. Orthlets are then formed as 
linear combinations of the shape functions and the added 
singularity functions. In three dimensions, we form shape 
functions as cartesian products of the ID shape functions, 
Sji(r) = S na: (x)S ny (y)S nz (z). In preliminary calcula- 
tions to represent the hydrogen atom ground state, we 
have found that adding a set of narrow gaussians, mul- 
tiplied by a windowing function similar to p(x), is very 
convenient and effective for representing the cusp singu- 
larity at the nucleus. Part of the convenience is that 3D 
gaussians are products of ID gaussians, like the 3D shape 
functions. These results will be presented elsewhere. 

There are several ways to change the basic lattice spac- 
ing in different regions of space. One simple approach is 
to let the finer and coarser grids overlap, so that com- 
pleteness is ensured, and then apply a GS procedure to 
the overlap region, automatically generating local func- 
tions which connect the two regions. However, the ability 
to generate orthlets with singularities means that chang- 
ing the lattice spacing is usually not necessary. 

Now, suppose one wishes to generate orthlets to repre- 
sent a set of functions {f a (r}}- For example, one might 
want to build an orthlet basis which is able to represent a 
standard radial basis set from quantum chemistry, since 
these basis sets are known to represent Hartree Fock or- 
bitals well. Note that the f a need not be orthogonal, but 
the orthlet basis generated from them is. The orthlet 
basis will automatically have additional degrees of free- 
dom allowing improved treatment of correlations. For 
simplicity, we will let Sft(r) represent both shape func- 
tions and any extra singularity basis functions. To adapt 
the basis to represent the f a , we apply the procedure in 
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the density matrix renormalization group for targetting 
more than one state |(|. For each lattice site j, we find 
the coefficients 

c a n = (Sl(m a (r))- (15) 
We now form the positive semidefinite density matrix 



n'n / j 



(16) 



Note that positive weighting factors a a can also be in- 
cluded in the sum if some of the functions are considered 
more important than others. The eigenvectors of p? , v~, 
define orthlet functions v^(r) — Y]^ v~Si(r) which opti- 
mally represent the functions {/°(r)} in the basis for site 
j. Each density matrix eigenvalue wp gives the weight 
associated with that orthlet in representing the f a . By 
choosing a cut off weight, and retaining all orthlets with 
weight greater than the cut off, one obtains a systemati- 
cally improvable basis for the f a . In particular, one can 
show that 



(17) 



Thus, if a density matrix eigenvalue wp is very small, 
then none of the f a have significant overlap with the 
corresponding function v@ , and v@ need not be included 
in the basis. The basis set becomes exact if the number of 
orthlets kept per site is equal to the number of functions 
f a . If only one function is in the set f a , then there is 
only one nonzero eigenvalue and the orthlet is simply the 
normalized projected function P t f a - 

As an example of this procedure, we generate a gen- 
eral set of orthlets to use as a basis set for the tails of 
wavefunctions. The ordinary shape functions are too lo- 
calized to represent an exponentially decaying tail effi- 
ciently. The orthlets we generate here extend quite far 
in one direction, and so are very useful to use as replace- 
ments for the shape functions for the edge sites of the 
lattice, allowing fewer sites to be used. The functions 
are constructed by using as f a a set of 13 gaussians with 
widths ranging from 1 to 4, each centered at the ori- 
gin, constructing and diagonalizaing the density matrix. 
The basis for an orthlet includes more than one site here: 
we include all 78 shape functions on sites 3 to 15. The 
resulting orthlets are able to represent any linear com- 
bination of the 13 gaussians with good accuracy. The 
three most important functions, along with the density 
matrix eigenvalues uup, are shown in Fig. 3. One can see 
the the density matrix eigenvalues decay very rapidly, so 
that only a few of these orthlets are needed. 

In electronic structure calculations, one would use a 
lattice spacing appropriate for the valence electrons, say 
0.3-1.0 angstroms. Cores would be treated using orthlets 
derived from localized cusp functions, which would be 
tied to each nucleus. 
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FIG. 3. Orthlets which can be used to represent tails of 
wavefunctions. The orthlets are adapted to a set of 13 gaus- 
sians centered at with a variety of widths. The density 
matrix eigenvalue w of each orthlet is given. The arrow in- 
dicates the location of the first site (3) included in the basis 
for the orthlets. The functions are orthogonal to all shape 
functions to the left of that site. 

Lattice sites near cores would have a dozen or more 
orthlets; in other areas, we expect only a few might 
be needed, making perhaps 100-300 basis functions per 
atom. Although this is an order of magnitude more than 
with the radial functions used in quantum chemistry, it 
is perhaps an order of magnitude less than a wavelet ba- 
sis, and the orthlets would be orthogonal and compact. 
The orthlets appear to be very convenient for the de- 
velopment of o(N) algorithms for density functional cal- 
culations; for example, sparse matrix methods coupled 
with the multigrid algorithm might be used to solve the 
Poisson equation. Orthlets also appear well adapted to 
multipole representations of the electron-electron inter- 
action, which can also be used in o(N) algorithms. 

We acknowledge support from the NSF under Grant 
No. DMR-9870930. 
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